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This paper presents a free-form deformation technique suitable for aerodynamic shape 
optimization. Because the proposed technique is independent of grid topology, we can treat 
structured and unstructured computational fluid dynamics grids in the same manner. The 
proposed technique is an alternative shape parameterization technique to a trivariate 
volume technique. It retains the flexibility and freedom of trivariate volumes for CFD shape 
optimization, but it uses a bivariate surface representation. This reduces the number of 
design variables by an order of magnitude, and it provides much better control for surface 
shape changes. The proposed technique is simple, compact, and efficient. The analytical 
sensitivity derivatives are independent of the design variables and are easily computed for 
use in a gradient-based optimization. The paper includes the complete formulation and 
aerodynamics shape optimization results. 
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Nomenclature 

B-spline basis functions 

product of B-spline basis functions 

maximum number of points in u-direction 

maximum number of points in v-direction 

number of surface points 

control points 

weights 

baseline grid 

grid point number 

coordinates (x, y, z) 

design variable vector 

surface bivariate coordinates 


I. Introduction 

Over the past several decades, aerodynamic shape optimization (ASO) has been successfully applied for two- 
dimensional and simple three-dimensional configurations. Four distinct steps are involved in ASO: 1) geometry 
parameterization; 2) surface grid generation; 3) volume grid generation, regeneration, or deformation; and 4) 
computational fluid dynamics (CFD) function and sensitivity analyses. Automatic grid generation and grid 
sensitivity analysis are two important issues for ASO. CFD grid generation for a full airplane model is time- 
consuming and costly. Detailed CFD grids based on a computer-aided design (CAD) model can take several weeks 
to develop and validate. Gradient-based optimization requires accurate sensitivity analysis for all four steps of ASO. 
Steps 2-4 can be automated for general configuration with unstructured grid technology and for specific 
configurations with structured grid technology. The level of effort for the geometric parameterization depends on the 
complexity of the vehicle configuration. This step is a labor-intensive process, and it is difficult to automate. For 
more details, readers are referred to an overview paper by this author, Ref. 1, on geometric modeling and grid 
generation for design and optimization. The purpose of this paper is to present a technique to simplify step 1 (shape 
parameterization) of ASO. 
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In general, shape parameterization techniques can be divided into eight categories 2 : basis vector, domain 
element, discrete, analytical, free-form deformation (FFD), partial differential equation (PDE), polynomial and 
spline, and computer-aided design (CAD). Among these categories, FFD and CAD appear to have the following 
desirable characteristics: 1) efficiency, 2) compactness, and 3) suitability for complex configurations. In addition, 
CAD-based parameterization techniques can achieve large geometry changes. This capability makes CAD-based 
techniques suitable for the conceptual design efforts. Another advantage of CAD-based optimization is the 
availability of a comprehensive set of geometric functionalities provided by commercial CAD systems. However, 
existing CAD systems lack tools for the automation of geometry abstraction and analytical sensitivity derivatives 
calculations. 

Commercial CAD systems use the feature -based solid modeling (FBSM) concept to create parametric design 
models. To parameterize an existing complex model is still a challenging task with today's CAD systems, and the 
models created are not always good enough for automatic grid generation tools. Even though the use of parametric 
modeling in design would make FBSM tools ideal for optimization, existing FBSM tools are not capable of 
calculating sensitivity derivatives analytically. Computer codes for commercial CAD systems are very large; to 
differentiate an entire system with automatic differentiation tools may not be a trivial task. Therefore, calculation of 
analytical sensitivity derivatives of geometry with respect to design variables could prove to be difficult within a 
commercial CAD environment. For some limited cases, the analytical shape sensitivity derivatives can be calculated 
based on a CAD model; however, this method will not work under all circumstances. One difficulty is that, for 
perturbation of some parameters, the topology of the CAD part may be changed. Another way to calculate the 
sensitivity derivatives is to use finite differences, as long as the perturbed geometry has the same topology as the 
unperturbed geometry. Both methods — the analytical and finite-difference approximations — have their difficulties 
and limitations for CAD-based parameterization. 

Free-form deformation is a subset of the soft object animation (SOA) 2 algorithms used in computer graphics for 
moiphing images and deforming models. The FFD technique can achieve only small to medium geometry changes. 
Because the FFD parameterizes the grid, it eliminates the need for geometry abstraction and surface grid generation. 
Unlike the CAD-based technique, FFD does not require one to reverse engineer the original design parameters. The 
FFD technique is suitable for the preliminary design phase, where geometry changes are rather small. Another 
advantage of parameterizing the grid is that the grid topology stays fixed throughout optimization; hence, the grids 
can be regenerated (deformed) automatically. The disadvantage is that the large shape changes could produce 
unacceptable grids due to the fixed topology. 

Our proposed technique uses nonuniform rational B-spline (NURBS) to model the geometry perturbations 
rather than the geometry itself. The parameterization of the shape perturbations is used in conjunction with the FFD. 
An appendix is included to provide a brief overview of NURBS representation and terminologies. 

II. Shape Perturbations 

During the preliminary design phase, numerical optimization starts with an existing design, and the goal is to 
redesign and improve performance by using numerical optimization. The geometric changes (perturbations) between 
initial and optimized shape are very small 2 , but the difference in the performance can be substantial. An effective 
way to reduce the number of shape design variables is to parameterize the shape perturbations instead of 
parameterizing the shape itself. Throughout the optimization cycles, the analysis grid can then be updated as 

r n( V ) = r !i + Ar n( V ) (!) 

where r n is the design grid point n, v is design variable vector, r* is the baseline grid, and A r n (v) is the shape 
design perturbation for grid point n. Far fewer design variables are required to parameterize the shape perturbations, 
A r (v) , than the baseline shape itself. Because the parameterization is independent of grid topology, we can use 
this technique for structured and unstructured CFD. 

III. Trivariate Volume-Based Deformation Techniques 

The use of a tri variate volume -based deformation technique for shape optimization has been discussed in Ref. 3. 
The FFD algorithm described by Sederberg and Parry 4 is ideal for deforming the polygonal models (e.g., CFD 
grids). Like other SOA algorithms, this algorithm maintains the polygon connectivity, and the deformation is 
applied only to the vertices of the model. The FFD process is analogous to embedding the grid inside a block of 
clear, flexible plastic (deformation object) so that, as the plastic is deformed, the grid is deformed as well. 
Deformation of complex shapes may require several deformation objects. The shapes of these deformation objects 
are not arbitrary. In fact, the shapes are trivariate volumes and could range from a parallelepiped to a general 
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NURBS volume as shown in Fig. la. The grid is deformed by perturbing the vertices that control the shape of the 
deformation block (e.g., corners of the parallelepiped). Sederberg and Parry used a Bezier volume 4 . Coquillart at 
I N R I A extended the Bezier parallelepiped to a nonparallelepiped cubic Bezier volume 5 . Lamousin and 
Waggenspack have further generalized this idea to NURBS volume 6 (Fig. lb), and they used multiple blocks to 
model complex shapes. This technique has been used for design and optimization by Yeh and Vance 7 , Perry and 
Balling 8 , and Ronzheimer 9 . For trivariate volume blocks, coordinates controlling the deformation are related through 
the mapping coordinates (£, T|, t). These coordinates are used in both forward and backward mapping. For more 
details, readers are referred to the paper by this author. Ref. 2. 

IV. Proposed Technique 

The trivariate volume operates on the whole space, regardless of the representation of the deformed objects 
embedded in the space (Fig. 2). Aerodynamic shape optimization is concerned with modifying the aerodynamic 
surfaces, not the volumes. Because the trivariate representation operates on the whole space, the design variables 
may have no physical significance for the design engineers. This drawback makes it difficult to select an effective 
and compact set of design variables. 

If we shrink the trivariate volume along one coordinate direction (e.g., Q, it will be reduced to a bivariate FFD 
surface (Fig. 3). These surfaces can be laid (pasted) onto outer mold line (OML) surfaces. Then as these bivariate 
FFD surfaces move, the underlying grid points move as well. It is much more effective to modify the OML surfaces 
by these bivariate FFD surfaces than by trivariate volumes. The proposed technique retains the flexibility and 
freedom of trivariate volumes for CFD shape optimization, but it uses a bivariate representation. This reduces the 
number of design variables by an order of magnitude, and it provides a much better control for surface shape 
changes. Because the NURBS formulation is the most general free-form surface representation, it is used for the 
bivariate surfaces in this paper. 

The implementation of the proposed algorithm is divided into the following steps: 

1. Identify regions of OML that need to be designed or adjusted (e.g., fillet, inboard leading edge). These 
regions are identified by pasting a set of NURBS surfaces on the outer mold line (OML). We will refer to 
these surfaces as marking surfaces. These surfaces are used to mark and identify surface grid points that 
belong to individual design regions. 

2. Determine the bivariate coordinates (q.q) of each grid point by projecting each one onto the marking 
surfaces. 

3. Create design surfaces corresponding to each marking surface defined in step 1. 

Step 1: Identify Design Regions 

CAD systems are frequently used during the preliminary design phase to model and represent the OML 
geometry. These representations are either bivariate functions (e.g., NURBS, Bezier) or implicit functions 
(cylinders, spheres), which can be converted to bivariate functions. Figure 4 shows a CAD representation of a 
business jet model. The solid lines are the boundaries of CAD surfaces. Any bivariate surface can be used to 
represent the marking surface. We have used a NURBS representation, but the algorithm presented here is equally 
applicable for any bivariate surface representation. Figure 4 shows a marking surface for a business jet. Two 
important requirements are needed for creating marking surfaces: 1) they cannot overlap each other, but they can 
overlap CAD surfaces, and 2) they must be laid very close to OML surfaces within some tolerance. The number of 
control points and order of marking surfaces play no role in the parameterization process. However, the knot vector 
distribution could play a role in the quality of shape deformation. 

Step 2: Project Grid Points on Marking Surfaces 

Each baseline surface grid point is projected to all marking surfaces created in step 1. Then, each grid point is 
linked to the closest marking surface. The algorithm for grid point projection can be found in Ref. 10. After this 
step, each grid point in the design region is assigned a design surface (closest surface) and corresponding bivariate 
coordinates (g ,Tj ). 

Step 3: Create NURBS Design Surfaces 

The Ar(v) in Eq. (1) is represented by a NURBS (see appendix for details on NURBS representation): 
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(2) 


M 


A r„ (v) = A r(v, £ , T ) n ) = X C », (£, » V n K 

m 


The coefficients C„, are products of B-spline basis functions. The only requirement for this design surface is that its 
knot vectors must span over bivariate coordinates r|) of all grid points belonging to individual design surfaces. 
We are free to pick the NURBS order and number of control points. The latter allows us to easily change the number 
design variable regardless of complexity of the baseline surfaces. As the control points of NURBS design surfaces 
move, the grid points belonging to the surface will move as well. The movement is based on an inverse mapping 
between the grid point bivariate coordinates and the NURBS surface. The bivariate coordinates of the grid points 
with respect to the marking surfaces are kept fixed through the optimization cycles. The movement may be 
unconstrained, or it can be constrained to move normally to the OML. The latter will reduce the number of control 
points (design variables) by a factor of three. All NURBS control points are initially set to zero, so the first 
computation gives the baseline grid. 

The analytical sensitivity derivatives of grid points with respect to the design variables are determined from 
Eq.(2). 


j^Cv) 

dv ,n 


d(Ar„(v)) _ 
dv 

m 


C,„ (^7i ’ Vn ) 


(3) 


The grid point sensitivity is independent of the design variables. Thus, we needed to calculate the sensitivity only 
once, at the beginning of the optimization. Then grid points can be updated as 


M 


(4) 


The C„, coefficients are precomputed, and they are used for sensitivity analysis as in Eq. (3) and for updating the 
grid point locations, as Eq. (4). 


Results 

Results are presented for a generic business jet and a morphing aircraft vehicle. Fig. 5 shows the deformed fillet 
for the business jet shown in Fig. 4. The marking surface was represented by a NURBS with 24 control points. The 
control points on the boundary were not allowed to move. The deformation was based on Eq.(4). 

The second example is for the aerodynamic shape optimization of a morphing aircraft vehicle. 11 The overall 
vehicle optimization requires multiple disciplines (e.g., structures, control) and constraints to provide useful results 
for design. The purpose of these optimization examples is only to demonstrate the shape parameterization technique 
presented in this paper. Half of the vehicle was used for the aerodynamic optimization (Fig. 6). The FUN3D 12 CFD 
program was used for sensitivity analysis. This program is an unstructured Euler and Reynolds-averaged Navier- 
Stokes solvers, which computes the sensitivity derivatives efficiently and accurately. FUN3D is a node-based, 
implicit, upwind, finite volume solver on mixed-element unstructured grids across the speed range. The discrete 
adjoint sensitivity derivatives are hand coded. The Euler mesh was generated for a Mach 0.8 and angle of attack of 
2.0. The mesh had 238362 nodes (1401823 tetrahedral cells). Fig. 7 shows the surface grid. Fig. 8 shows the 
marking NURBS surfaces used for the shape parameterization. There are 150 control points (design variables) 
available for optimization, but 54 design variables were active for shape optimization. Fig. 9 shows the surface grid 
points that are allowed to move. Two optimization calculations were performed on a 32-CPU Pentium IV cluster. 
The first optimization run was a lift-constrained drag minimization, resulting in a 23 percent reduction in drag (Fig. 
10) while maintaining constant lift. The second case is an unconstrained L/D maximization, resulting in a 32 percent 
increase in L/D (Fig. 11). It took less than a day to prepare the CAD model, generate the Euler grid, set up the 
parameterization model, and set up CFD runs and optimization. It took another day to complete the calculations. 

Conclusion 

This paper presented a free-form deformation technique suitable for aerodynamic shape optimization. The 
present technique reduced the number of design variables by an order of magnitude. The analytical sensitivity 
derivatives were independent of the design variables and were easily computed for use in a gradient-based 
optimization. The aerodynamic shape optimization results demonstrated the effectiveness of the present formulation. 
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Appendix 

This section contains a brief overview of the NURBS, and readers should consult Ref. 13 and 14 for a detailed 
discussion. Before 1 describe a NURBS surface representation, I will present a brief description of NURBS curves. 
A NURBS curve, r(g) , can be represented as 

Uf) = (A.l) 

1=1 

The univariate coordinate, q, is bounded by ^ min < t, < £ max . The V, are the control points (forming a control 
polygon), and W, are the weights. The B ip are the p - th degree B-spline basis functions defined on the nonperiodic 
and nonuniform knot vector (£) 

4 G [4mm = ^ ^ ^ 4 = 4max ] ( A - 2 ) 

where k is the number of knots. This completes the mapping between the one-dimensional univariate coordinate, 
and the three-dimensional Euclidean space, r. A NURBS curve has five important properties: 
o It is invariant under linear transformation. 

o A NURBS curve of order p, having no multiple interior knots, is p - 2 differentiable, 
o The approximation is local in nature. 

o A NURBS curve is contained in the convex hull of its control points, 
o The NURBS approximation is variation diminishing. 

To evaluate, the three-dimensional curve NURBS is commonly represented in homogeneous form as 

r(x,y,x) <=> r w (wx,wy,wz,w) (A.3) 

So, the NURBS curve can conveniently be defined as a perspective map of its nonrational counterpart in four- 
dimensional space as 

< a - 4 > 

M 

where V ]" is defined as V-" = {x i W j ,y j W i , z.w t , W ; } . The basis functions can be efficiently computed using 
DeBoor’s algorithm 13 . 

A NURBS surface is a parametric surface and is defined as a bivariate function as 

'is. n) = yyj (a. 5) 

i j 

where the Vy are the locations of NURBS control points, the W ;j coefficients are the weights, £, and T] are the 
bivariate coordinates, / and J are the numbers of control points in ^ and !] directions, and B n (^) and 
B . q (l]) are the B-spline basis functions of degrees p and q respectively. 

Similar to a NURBS curve, the NURBS surface can conveniently be defined as a perspective map of its 
nonrational counterpart in four-dimensional space as 

= (A.6) 

i j 

This equation can be written in a compact form by collapsing the two summations into a single summation as 

M 

r (£;,Jl) = Y J C m (^B)V*\ where m = i + I*(j- 1), M = I*J, and 

m 

cM,n)=c t y,n) = s,,(#>s p ,w <a.?> 
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Figure 1: a) Parallelepiped volume, b) NURBS volume. 



a) 



b) 

Figure 2: a) Trivariate volume deformation, b) Deformed shape. 
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a) 



b) 


Figure 3: a) Bivariate surface deformation, b) Deformed cross section. 



Figure 4: Parameterizing wing-fuselage fillet of a bushiness jet. 
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Figure 5: Free-form deformation of wing-fuselage fillet. 



Figure 6: Morphing vehicle. 
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Figure 7: Surface grid. 



Figure 8: NURBS marking surfaces. 



Figure 9: Surface grid points identified for deformation. 
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Figure 11: Unconstrained L/D maximization. 
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